type triangle_tree =
	| Null
	| Node of (triangle_tree*Geometry.trianglei*triangle_tree)

(* to_lin
 * Convert matrix2 coordinates to linear coordinates.
 * Inputs:
       * w as int : with of the matrix.
       * x as int : x coordinate.
       * y as int : y coordinate.
 *)
let to_lin w x y =
    y * w + x

(* calc_z
 * Caculates altitude of a given point.
 * Inputs:
       * hinit as float -> base altitude.
       * href as float -> reference altitude.
       * dist as float -> dist to nearest border point.
       * tau as float -> magic param.
 *)
let sq2 = 1.
let calc_z (hinit,href) distance tau scale =
  let dist = distance *.scale
  in
    if dist>=5.*.tau then
      hinit
    else if dist<=sq2 then
      href
    else
      hinit -. (hinit-.href)*.(exp (-.dist/.tau))

(*
     let dh = hinit -. href in
      hinit -. dh*.(exp (-.dist/.tau))
 *)

(* get_implicit_triangle
 * Calculates implicit triangle acording to the triangle's
 * implicit grid-coordinate.
 * Inputs:
      * x as int = x coordinate of a square of the grid.
                   Points to the n'th cell where n = x div 2.
                   An even x points to the left sub-triangle of
                   the square. An odd x represents the right one.
      * y as int = y coordinate of a square of the grid.
      * cw as float = width of a cell
      * ch as float = height of a cell
 * Outputs:
      * val as trianglei = one of the two sub-triangle of the cell
 *)
let get_implicit_triangle x y lf=
   let xdiv2 = x/2 in
      {Geometry.isa=lf xdiv2 y;
       Geometry.isb=lf (xdiv2+1) (y + ((x+1) mod 2));
       Geometry.isc=lf (xdiv2 + (x mod 2)) (y+1)}


(* write_tr
 * Writes a triangle to a previously opened obj file.
 * Inputs:
      * i as int ref = reference to a count of written triangles
      * t as file descriptor = opened obj's file handle. The file
                               must be opened with write-access.
 *)
let write_face file t =
   ObjExporter.write_face file t
let write_point file p map scale =
  map.(int_of_float p.Geometry.x).(int_of_float p.Geometry.y)
  <- Some(p.Geometry.z);
  ObjExporter.write_point file p scale

(* map_implicit_triangles
 * Maps two exclusives functions on each implicit triangle of an
 * implicit grid.
 * Inputs:
       * b as int -> int -> int ref -> bool =
            first function to map to each intersection point of
            the grid.
       * f as int ref -> trianglei -> unit =
            function called iff b returns false. The triangle
            pushed to f is the current grid-cell implicit triangle.
       * p as int -> size of an edge of a cell.
       * nx as int -> number of cell on the x axis.
       * ny as int -> numbre of cell on the y axis.
 *)
let map_implicit_triangles lf b f nx ny =
   for j=0 to ny-1 do
      for i=0 to nx-1 do
         if not (b i j) then
            begin
	    f (get_implicit_triangle i j lf)
            end
      done
   done

(* leaf_map
 * Applies a given function to each leaf'value of a triangle_tree.
 * Inputs:
	* f as triangle -> unit : function to apply
	* _ as triangle_tree : tree to explore
 *)
let rec leaf_map f = function
   Null -> ()
   | Node(fg,v,fd) ->
	match fg with
	   Null -> f v
	   | Node(_) ->
                begin
                leaf_map f fg;
                leaf_map f fd;
                end

(* insert_childs_f
 * Inserts a couple of child on the leaf selected
 * by predicate predicate.
 * Inputs:
        * predicate as trianglei -> bool : predicate which
                                    determine which child
                                    have to be selected.
        * f as trianglei -> (triangle_tree, trianglei, triangle_tree) :
                            function splitting a node into two childs.
        * val as triangle_tree : tree to explore.
  *)
let rec insert_childs_f predicate f =function
   Null -> Null
   | Node(Null,v,Null) -> f v
   | Node(fg,v,fd) ->
         begin
         match fg with
            Node(_,vfg,_) ->
               if (predicate vfg) then
                  Node(((insert_childs_f predicate f fg),v,fd))
               else
                  Node((fg,v,(insert_childs_f predicate f fd)))
           | Null -> Null
          end

(* tree_map
 * Applies a given function to each node's value of a triangle_tree.
 * Inputs:
	* f as triangle -> unit : function to apply
	* _ as triangle_tree : tree to explore
 *)
let rec tree_map f = function
   Null -> ()
   | Node(fg,v,fd) ->
         if f v then
            tree_map f fg
         else
            tree_map f fd

(* write_infected_triangles
 * Write an infected triangle tree's leafs to the obj.
 * Inputs:
       file as output chanel : opened obj file.
       leaf_t as triangle_tree array array : trangle tree array to explore.
       x as int : x coordinate of the tree in leaf_t
       y as int : y coordinate of the tree in leaf_t
 *)
let write_infected_triangles file leaf_t x y =
   (* test if the triangle is infected *)
   if leaf_t.(x).(y)<>Null then
      begin
         leaf_map (write_face file) leaf_t.(x).(y) ;
         true;
      end
   else (* not infected: no triangle to write *)
      false

(* kill_0disk
 * Sets a rectangle centered on a given points to -1.
 *)
let kill_0disk l x y r w h=
   for i = x - r to x + r do
      for j = y - r to y + r do
         if i < w && j < h && i>=0 && j>=0 then
           l.(i).(j) <- -1.
      done
   done

(* pt_in_rect
 * Calculates splited edge given by the intersection between a
 * point and a couple of triangles (i.e a rectangle)
 * Inputs:
       pt as point array -> points
       p as float -> grid precision
       i as int -> x : coordinate of the splitting point
       j as int -> y : coordinate of the splitting point
       to_lin_fun as int -> int -> int :
              matrix-coordinates to linear coordinates function.
       pt_iterator as int -> point : iterator of a points matrix.
 *)
let pt_in_rect pt p i j to_lin_fun pt_iterator=
      let ip = i/p*2 and
          jp = j/p in
      let t0 = get_implicit_triangle ip jp to_lin_fun in
      match Geometry.pti_in_tiborders pt_iterator t0 pt with(*triangle:|\*)
         | 0 ->
	    let t1=  get_implicit_triangle (ip+1) jp to_lin_fun in
            (match Geometry.pti_in_tiborders pt_iterator t1 pt with(* \| *)
                   | 0 -> ((-1,-1),(-1,-1),-1,-1)
                   | 1 -> (* _ *) ((ip,jp-1),(ip+1,jp),2,1)
                   | 2 -> (* | *) ((ip+2,jp),(ip+1,jp),3,2)
                   | _ -> (* \ *) ((ip,jp),(ip+1,jp),1,3)
              )
         | 1 -> (* \ *) ((ip+1,jp),(ip,jp),3,1)
         | 2 -> (* _ *) ((ip+1,jp+1),(ip,jp),1,2)
         | _ -> (* | *) ((ip-1,jp),(ip,jp),2,3)


(* infect_triangle
 * Infects a given triangle, splitting it in two
 * and inserting it in the correspouding
 * binary tree.
 *)
let infect_triangle (x1,y1) grid_points iedge
                    nx  i j m dm
                    infected_l =
   if infected_l.(x1).(y1) = Null then (* not yet infected *)
      let t = get_implicit_triangle x1 y1 (to_lin (nx+1)) in
      let (t1,t2) =
         Geometry.trianglei_split_segi_n t
                                      ((DynArray.size grid_points)-1)
                                      iedge in
         infected_l.(x1).(y1)<- Node(
                                  (Node(Null,t1,Null)
                                  ,{Geometry.isa=t.Geometry.isa;
                                    Geometry.isb=t.Geometry.isb;
                                    Geometry.isc=t.Geometry.isc;}
                                  ,Node(Null,t2,Null)))
   else (* already infected! Infect final node *)
      begin
          let predicate = (function t ->
                             (Geometry.pti_in_tiborders
                             (DynArray.get grid_points)
                              t
                              {Geometry.x=float_of_int i;
                               Geometry.y=float_of_int j;
                               Geometry.z=0.}
                              )) in
           let funct triangl=
                    begin
                       let (t1,t2) = Geometry.trianglei_split_segi_n
                             triangl ((DynArray.size grid_points)-1)
                                     (predicate triangl);
                       in
                       Node((Node((Null,t1,Null))),
                            triangl,
                            (Node((Null,t2,Null))))
                    end
             in
             infected_l.(x1).(y1)<- insert_childs_f
                                       (function t -> (predicate t) <> 0)
                                       funct
                                       infected_l.(x1).(y1);
   end

(* get_precision
 * Read precision file and returns datas.
 * Inputs:
       path as string : File path.
 *)
let get_precision path =
   try
   let f = open_in path in
   let hardautosmooth = (int_of_string (input_line f)) = 1 in
   let autosmooth = (int_of_string (input_line f)) = 1 in
   let tau =  float_of_string (input_line f) in
   let p = int_of_string (input_line f) in
       close_in f;
    (hardautosmooth, autosmooth, tau, p)
   with _ ->
      let f = open_out path in
         Printf.fprintf f "%s\n" "0";
         Printf.fprintf f "%s\n" "1";
         Printf.fprintf f "%s\n" "2.";
         Printf.fprintf f "%s\n" "5";
         close_out f;
         (false,true,2.,5)

(* write_triangles
 * Write the full implicite-triangle matrix to a graphic
 * file given by the graphics engine. Alse takes care of
 * infected triangles.
 * Inputs:
	* m as float array: distance map.
        * dm as float * float array: grid of
                                    (reference altitude * initial altetude)
 *)
exception TesselatorExp
let write_triangles_with_infected path m dma w h tau p scale =
  try
    let dm = Array.init (Array.length dma) (fun i -> Array.copy dma.(i))
    in
   let map = Array.init h (fun _ -> Array.make w None) in
   let f = ObjExporter.openObj path in
   let (hardautosmooth, autosmooth) = (false,true) in
   let modmod = if (w mod p = 0) then 1 else 0 in
   let nx = w/p - modmod and
       ny = h/p - modmod in
   let infected_l = Array.init (nx*2) (fun _ -> (Array.make ny Null)) in
   let grid_points = DynArray.make ((nx+1)*(ny+1))
                                   {Geometry.x=0.;
                                    Geometry.y=0.;
                                    Geometry.z=0.} in
       begin
          (* Construct grid_points grid explicitly (with altitude)
           * and kill infected points in grid intersections *)
          for j=0 to ny do
             for i=0 to nx do
                   DynArray.set grid_points
                      ({Geometry.x=float_of_int (i*p);
                        Geometry.y=float_of_int (j*p);
                        Geometry.z=(calc_z m.(i*p).(j*p)
                                    dm.(i*p).(j*p)
                                    tau scale)})
                        (to_lin (nx+1) i j);
  	        if dm.(i*p).(j*p) = 0. then
		   kill_0disk dm (i*p) (j*p) 1 w h;
             done
          done;
          if autosmooth then
             begin
             (* Mark infected triangles; kill disk; and build b-tree *)
             for j=0 to h-1 do
                for i=0 to w-1 do
                  if (hardautosmooth &&
                      dm.(i).(j) >= 0. &&
                      dm.(i).(j) <= 5.*.tau)
                     || dm.(i).(j) = 0. then
                  begin
                     let ((x1,y1),(x2,y2),ia,ib)=pt_in_rect
                                           {Geometry.x=float_of_int i;
                                            Geometry.y=float_of_int j;
                                            Geometry.z=0.}
                                           p i j (to_lin (nx+1))
                                          (DynArray.get grid_points) in

                     if ib<>(-1) then
                     begin

                     let cnd1 = x1>=0 && x1<nx*2 && y1>=0 && y1<ny in
                     let cnd2 =  x2>=0 && x2<nx*2 && y2>=0 && y2<ny in
                     if (cnd1 || cnd2) then
                     begin
                          DynArray.add grid_points
                                   {Geometry.x=float_of_int i;
                                    Geometry.y=float_of_int j;
                                    Geometry.z=calc_z m.(i).(j)
                                                      dm.(i).(j)
                                                      tau scale};
                      end;
                      if cnd1 then
                        infect_triangle (x1,y1) grid_points
                                        ia nx i j m dm
                                        infected_l;

                     if cnd2  then
                        infect_triangle (x2,y2) grid_points
                                        ib nx i j m dm
                                        infected_l;
                     if (cnd1 || cnd2) then
                        kill_0disk dm i j 1 w h
                    end
                  end
                done;
             done;
         end;
         (* Write points *)
         let s = DynArray.size grid_points in
             for i=0 to s-1 do
                write_point f (DynArray.get grid_points i) map scale
             done;

         (* Write triangles (faces) *)
         map_implicit_triangles (to_lin (nx+1))
                                (write_infected_triangles f infected_l)
                                (write_face f)
                                (nx*2) ny;

         (* Close file *)
         ObjExporter.closeObj f;
         map
      end
  with
      _-> raise TesselatorExp

